#ifndef SPECIES3D_H
#define SPECIES3D_H

#include <vector>
#include <string>

#include "../Species/Species.h"

using namespace std;

class Species3D : public Species
{
public:
    //constructor
    Species3D(PicParams* params_, ParallelMPI* pmpi_, int ispec);

    //destructor
    ~Species3D();

    //calculate the particle dynamics (interpolation, pusher, projection)
    //only caculate the number density, no electric currents
    virtual void move(PhysicalField* fields);

    //project particles to grid points, to calculate charge density
    virtual void project(PhysicalField* fields);

    //interpolate fields at particle position, and move particle
    void push(CellParticles<3>* cell_particles, const LocalFields E_cell_fields[2][2][2], const LocalFields B_cell_fields[2][2][2]);

    //decide if particles hit boundary, cross MPI boundary, to move to another cell, or still stay in the cell
    virtual void sort_particles();

    //After MPI exchange particles, delete useless particles
    virtual void correct_particles();

    inline void interpolate(const LocalFields cell_fields[2][2][2], LocalFields& local_field, const double position[3])
    {
        //Normalized particle position
        double xpn = position[0] * params->cell_length_inv[0];
        double ypn = position[1] * params->cell_length_inv[1];
        double zpn = position[2] * params->cell_length_inv[2];

        //Indexes of the central nodes
        int ip_ = floor(xpn);
        int jp_ = floor(ypn);
        int kp_ = floor(zpn);


        //Declaration and calculation of the coefficient for interpolation
        double delta;
        double coeffxp[2], coeffyp[2], coeffzp[2];

        delta   = xpn - (double)ip_;
        coeffxp[0] = 1.0 - delta;
        coeffxp[1] = delta;


        delta   = ypn - (double)jp_;
        coeffyp[0] = 1.0 - delta;
        coeffyp[1] = delta;

        delta   = zpn - (double)kp_;
        coeffzp[0] = 1.0 - delta;
        coeffzp[1] = delta;

        local_field.x = 0.0;
        local_field.y = 0.0;
        local_field.z = 0.0;

        for (int iloc=0 ; iloc<2 ; iloc++)
        {
            for (int jloc=0 ; jloc<2 ; jloc++)
            {
                for (int kloc=0 ; kloc<2 ; kloc++)
                {
                    local_field.x += coeffxp[iloc] * coeffyp[jloc] * coeffzp[kloc] * cell_fields[iloc][jloc][kloc].x;
                    local_field.y += coeffxp[iloc] * coeffyp[jloc] * coeffzp[kloc] * cell_fields[iloc][jloc][kloc].y;
                    local_field.z += coeffxp[iloc] * coeffyp[jloc] * coeffzp[kloc] * cell_fields[iloc][jloc][kloc].z;
                }

            }
        }
    };

    virtual int get_number_of_particles()
    {
        int n = 0;
        for(int i = 0; i < particles_.size(); i++)
        {
            for(int j = 0; j < particles_[i].size(); j++)
            {
                for(int k = 0; k < particles_[i][j].size(); k++)
                {
                    n += particles_[i][j][k]->size();
                }
            }
        }
        return n;
    }

    void init_particles();

    vector< vector< vector< CellParticles<3>* > > > particles_;

    vector< Particle<3> > psi_particles_;

    vector< Particle<3> > exchange_particles_;

    PartBoundCond<3> part_bound_cond;

};



#endif
